library(rgdal)
library(plotly)
library(ggplot2)
library(dismo)
library(dplyr)
library(insol)
library(ggmap)
load("/home/rstudio/morph/data/test.rob")
Lat<-median(sites$lat) ## Get lat and long which are used in some functions such as day length
Lon<-median(sites$lon)
sites$mean_biomass[is.na(sites$mean_biomass)]<-0 ## Set not available to zero
map<-gmap(grat)
plot(map)
map2<-get_map(location=c(Lon,Lat),maptype="terrain",zoom=10)
map2<-ggmap(map2)
map2
FMakeTime<-function(year=2016,month=1,day=1,hr=1){
tm<-sprintf("%04d-%02d-%02d %02d:00:00",year,month,day,hr)
tm<-as.POSIXct(tm)
tm
}
FMakeDate<-function(ftm=tm)format(ftm,"%Y-%m-%d")
tm<-FMakeTime()
tm
## [1] "2016-01-01 01:00:00 UTC"
FMakeDate(tm)
## [1] "2016-01-01"
tm+60*60
## [1] "2016-01-01 02:00:00 UTC"
Use the insol package. This produces an object with sunrise and sunset times, so if the hour falls between them it is day.
FIsDay<-function(tm,Lat=55.32,Lon=-162.8)
{
hr<-as.numeric(format(tm, format='%H'))
day_len<-data.frame(daylength(Lat, Lon,JD(tm), tmz=-10))
isday<-ifelse(hr>day_len$sunrise & hr< day_len$sunset,"Day","Night")
isday}
tm<-FMakeTime(2016,1,1,10)
FIsDay(tm)
## [1] "Day"
tm<-FMakeTime(2016,1,1,7)
FIsDay(tm)
## [1] "Night"
This is the function taken from Humbolt paper
FBMR <- function(mass)4.59*mass^0.69
FBMR(1.4)
## [1] 5.7895
mass<-(1440:1576)/1000
plot(FBMR(mass)~mass)
FFat2Energy<-function(fat)34.3*fat # g fat to KJoules
FEnergy2Fat<-function(energy)energy/34.3 # KJoules to g fat
As I understand it from the spreadsheet this set of equation should calculate the metabolic rate in the same units ()
FBMR<-function(ftemperature=-10,fwindspeed=2,fmass=1500)
{
TBrant<-7.5
ftemperature[ftemperature>TBrant]<-TBrant
fwindspeed[fwindspeed<0.5]<-0.5
DeltaT<-TBrant-ftemperature
b<-0.0092*fmass^0.66*DeltaT^0.32
a<-4.15-b*sqrt(0.06)
a+b+sqrt(fwindspeed)
}
FBMR(10)
## [1] 5.564214
I haven’t got data for the whole of 2016.
clim<-subset(clim,as.numeric(format(clim$date,'%Y'))==2015)
temperature<-(clim$tmin+clim$tmax)/20
windspeed<-clim$avwind/10
plot(temperature,FBMR(temperature))
plot(clim$date,FBMR(temperature,windspeed))
biomass<-1:100
FConsumption<-function(fbiomass)100*0.01028*(1-exp(-0.105*fbiomass))*(1.0373*(1-exp(-0.0184*fbiomass))/60)
plot(FConsumption(biomass)~biomass)
Converting into energy assimilated
FEnergyAssim<-function(consumption, a=0.464, E= 16.8)consumption*a*E*1000
plot(FEnergyAssim(FConsumption(biomass))~biomass)
plot(FEnergyAssim(FConsumption(sites$mean_biomass))~sites$mean_biomass)
0.8.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1
Assume 20% over BMR. Add time argument in seconds. Returns energy used in KJoules.
BMR<-FBMR(temperature,windspeed,fmass=1700)
FRest<-function(fbmr=BMR,ftime=3600)
{1.2*fbmr*ftime/1000}
hist(FEnergy2Fat(FRest()))
Assume use twice BMR while feeding.
FFeed<-function(fbmr=BMR,ftime=3600,fbiomass=100){
EGain<-FEnergyAssim(FConsumption(fbiomass))*ftime/1000
EUse<-2*fbmr*ftime/1000
EGain-EUse
}
hist(FEnergy2Fat(FFeed()))
According to the Humbolt paper this is given by
FMaxDaily<-function(mass)1713*mass^0.72
FMaxDaily(1500)/1000
## [1] 331.5454
FFeed()[1]
## [1] 359.3228
So using these formulae the birds can get enough energy for a whole day from just one hour’s intensive feeding. This seems much too high. However this assumes that an hour feeding is completely dedicated to intensive feeding activity, which is unrealistic.
Assume velocity of 60km per second and overhead for take off for all flights equivalent to a minute’s flight time. Use 12 times BMR.
FFlightTime<-function(fspeed=60,fdistance=1000000)
{
speedms<-fspeed/3.6
ftime<-fdistance/speedms+60
}
FFly<-function(fbmr=BMR,ftime=FFlightTime()){
EUse<-12*fbmr*ftime/1000
EUse
}
hist(FEnergy2Fat(FFly()))
This function is the first in the set to receive two dataframes as arguments. The sites data frame is merged with the tides data frame after filtering the tides data frame to the rows representing the height of the tides in each bay. The proportion of the area in each patch lying between a minimum depth and a maximum height after changing the water level at each site accordingly is calculated, added to the sites data frame and the whole data frame is returned. So we need to assign the result to the sites data frame.
FSuitable<-function(fsites=sites,ftm=tm,ftides=tides,depth=-0.4,height=0){
current_tide<-subset(ftides,ftides$time==ftm)
d<-merge(fsites,current_tide)
tide<-d$ht
depth<-depth+tide
height<-height+tide
dd<-cbind(d$min,d$q10,d$q25,d$median,d$q75,d$q90,d$max,depth,height)
f<-function(x)
{
q<-c(0,10,25,50,75,90,100)
qs<-x[1:7]
depth<-x[8]
height<-x[9]
x2<-q[qs>=depth&qs<=height]
# The supress warnings are needed as the vector may be of zero
# length. This also leads to results of -inf instead of zero
suppressWarnings(x2<-max(x2,na.rm=TRUE)-min(x2,na.rm=TRUE))
if(is.na(x2))x2<-0
if(x2==-Inf)x2<-0
x2
}
fsites$psuitable<-apply(dd,1,f)
fsites
}
tm<-FMakeTime(2015,7,15)
dt<-FMakeDate(tm)
sites<-FSuitable(sites,tm)
str(sites)
## 'data.frame': 553 obs. of 21 variables:
## $ x : num -18134280 -18118072 -18114830 -18118072 -18129417 ...
## $ y : num 7424290 7435016 7435016 7422758 7424290 ...
## $ geom : num 0 0 0 0 0 0 0 0 0 0 ...
## $ area_m2 : num 805904 803687 803687 806221 805904 ...
## $ lon : num -163 -163 -163 -163 -163 ...
## $ lat : num 55.3 55.4 55.4 55.3 55.3 ...
## $ min : num -0.652 -7.435 -5.449 -0.75 -1.5 ...
## $ q10 : num 0.825 -4.164 -2.604 -0.75 -1.5 ...
## $ q25 : num 1.22 -1.876 -0.671 -0.585 -0.513 ...
## $ median : num 1.22 -0.715 -0.293 -0.105 0.553 ...
## $ mean : num 3.125 -1.445 -0.681 0.012 0.109 ...
## $ q75 : num 5.7275 -0.4179 -0.0661 0.5023 0.8234 ...
## $ q90 : num 7.22 -0.2411 0.0759 0.8652 0.8924 ...
## $ max : num 9.22 0.084 0.9 4.763 0.91 ...
## $ rid : int 251 106 120 285 264 497 452 496 455 209 ...
## $ psuitable : num 0 15 0 0 0 0 25 0 25 10 ...
## $ median_biomass : num 134 227 214 215 126 ...
## $ mean_biomass : num 131 213 201 183 131 ...
## $ median_shootlength: num 18 63.5 56 52 35 96 61 49 55 68 ...
## $ mean_shootlength : num 13.2 65.6 63.6 49.1 46 ...
## $ station : int 1 1 1 1 1 4 1 4 1 1 ...
The dataframes obtained from the data base include site data, distances between sites and climate data. However there are no birds. This function forms an initial birds data frame. Other derived characteristics can be added when the model runs.
FArriveBirds<-function(ftm=tm,nbirds= 10,fsites=sites,male_wt=1500,female_wt=1400,fat=300)
{
sex<-sample(c("M","F"),nbirds,replace=T)
lean_wt<-numeric(nbirds)
lean_wt[sex=="M"]<-male_wt
lean_wt[sex=="F"]<-female_wt
lean_wt<-lean_wt+rnorm(nbirds,0,sd=20)
fat<-rlnorm(nbirds,mean=log(fat),sd=0.2)
wt<-lean_wt+fat
energy_store<-FFat2Energy(fat)
## Add other properties here
##
rid<-sample(fsites$rid,nbirds,replace=TRUE) ## Place them at random
bid<-1:nbirds ## ID number
birds<-data.frame(bid,arrive_time=ftm,sex=sex,weight=wt,fat=fat,energy_store=energy_store,rid=rid)
birds
}
If we want to add birds after the first arrival time then we need to pass the birds dataframe that was created to the function in order to keep the ids sequential. The same function is used again to create more birds at the site.
FAddBirds<-function(ftm=tm,nbirds= 10,fsites=sites,fbirds=birds)
{
newbirds<-FArriveBirds(ftm,nbirds,fsites)
newbirds$bid<-newbirds$bid+max(fbirds$bid)
rbind(fbirds,newbirds)
}
birds<-FArriveBirds(tm,10,sites)
birds<-FAddBirds(tm,10000,sites,birds)
str(birds)
## 'data.frame': 10010 obs. of 7 variables:
## $ bid : int 1 2 3 4 5 6 7 8 9 10 ...
## $ arrive_time : POSIXct, format: "2015-07-15 01:00:00" "2015-07-15 01:00:00" ...
## $ sex : Factor w/ 2 levels "F","M": 2 2 1 1 1 1 1 1 2 2 ...
## $ weight : num 1808 1885 1661 1795 1715 ...
## $ fat : num 311 352 242 389 321 ...
## $ energy_store: num 10651 12090 8318 13350 11004 ...
## $ rid : int 351 307 470 429 355 489 223 334 368 8 ...
This is a tricky function to get right. The “dist”" data frame looks like this …
str(dist)
## 'data.frame': 305809 obs. of 3 variables:
## $ rid : int 251 251 251 251 251 251 251 251 251 251 ...
## $ rid2 : int 251 106 120 285 264 497 452 496 455 209 ...
## $ dist_m: num 0 11069 12651 9285 2773 ...
We have two rids, one representing the source site and the second representing the potential destinations. This contains all possible moves. To speed things up we only want to look at distances between sites with birds and only within a search distance.
FFilterMoves<-function(fdist=dist,fbirds=birds,search_distance=10000)
{
moves<-fdist[fdist$rid %in% unique(fbirds$rid),]
moves<-moves[moves$dist_m<search_distance,]
moves}
moves<-FFilterMoves()
head(moves)
## rid rid2 dist_m
## 1 251 251 0.0000
## 4 251 285 9284.9133
## 5 251 264 2772.8898
## 10 251 209 8717.6461
## 11 251 276 6527.9176
## 17 251 253 871.9961
We need to look at eelgrass biomass and the area that is suitable on the destination sites (psuitable). If we merge the sites data frame by default the matched column is rid, not rid2. I.e. the source rather than the destination. We should only look at sites with a suitability of greater than zero.
We can calculate the energy gained feeding and lost by moving using a fixed basal metabolic rate for comparative purposes to simplify a bit. There’s nothing to be gained by calculating the actual BMR on the day as all we really need at the end is a rank. The ranks are grouped by source rid using tapply.
FScoreMoves<-function(fmoves=moves,fsites=sites,flight_cost=TRUE)
{
fsites<-fsites[fsites$psuitable>0,]
fmoves<-merge(fmoves,fsites, by.x="rid2",by.y="rid")
fmoves$flight_time<-FFlightTime(fdistance=fmoves$dist_m)
fmoves$fly<-FFly(fbmr=9,ftime=fmoves$flight_time)
biomass<-fmoves$mean_biomass#*(fmoves$psuitable/100)
if (flight_cost==TRUE)
{fmoves$feed<-FFeed(fbmr=9,ftime=3600-fmoves$flight_time,fbiomass=biomass) ## Subtract the time neeed to fly there
fmoves$energy<-fmoves$feed -fmoves$fly}
if (flight_cost==FALSE)
{fmoves$feed<-FFeed(fbmr=9,ftime=3600,fbiomass=biomass) ## Subtract the time neeed to fly there
fmoves$energy<-fmoves$feed}
fmoves<-fmoves %>%
arrange(rid, -energy) %>%
group_by(rid) %>%
mutate(rank=row_number()) %>%
filter(rank==1)
fmoves<-fmoves[c("rid","rid2","dist_m","flight_time","feed","fly","energy","rank")]
}
moves<-FScoreMoves()
There will be some cases in which there are no suitable sites in range as they are all covered by the tide. In this case the birds stay put. Change the rid to the new value and add only the flight time and distance to the data frame.
FMoveBirds<-function(fbirds=birds,fmoves=moves)
{
bird_moves<-merge(fbirds,fmoves,all.x=TRUE)
bird_moves$rid2[is.na(bird_moves$rid2)]<-bird_moves$rid[is.na(bird_moves$rid2)]
fbirds$rid<-bird_moves$rid2
fbirds$distance<-bird_moves$dist_m
fbirds
}
moves<-FFilterMoves(search_distance=20000)
moves<-FScoreMoves()
birds<-FMoveBirds()
birds_sites<-merge(birds,sites)
map2+geom_point(data=birds_sites,aes(lon, lat),col="red")
for (i in 1:10)
{
tm<<-tm+(60*60*60)
sites<<-FSuitable(sites,tm)
moves<<-FFilterMoves(fbirds=birds)
moves<<-FScoreMoves()
birds<<-FMoveBirds(birds,moves)
birds_sites<-merge(birds,sites)
print(map2+geom_point(data=birds_sites,aes(lon, lat),col="red"))
}
save(list=ls(),file="/home/rstudio/morph/scripts/functions.rob")